function out = ResistanceABMFinal(sim,p)
% Models Nonviolent Resistance.

%%% Variables
% General Variables
a = .5; %Used in calculating income distribution
b = .5; %Used in calculating income distribution
k1 = 2.3; %Used in calculating probability of arrest
k2 = 0.025*3; % Used in calculating cost. I included the 3 so I could keep the same cost but drop the max jail term from 30 to 10 in order to speed up simulation.
if nargin < 2
    tic %Start time measurement for simulation
    % Simulation Variables
    sim.PauseTime = 0; % Seconds to pause each time step
    sim.PlotOn = 1; % 1 to make plots, 0 for no plots
    % Parameters
    p.MaxSteps = 200; % Maximum number of time steps
    p.LatticeX = 40; % Number of spaces in graph in X direction
    p.LatticeY = 40; % Number of spaces in graph in Y direction
    p.torus = 1; % Should lattice connect at sides? 1 Yes, 0 No, Currently, model is torus regardless
    p.PercentFillCitizens = 70; % Percent Lattice is filled with agent type
    p.PercentFillActivists = .8; % Percent Lattice is filled with agent type
    p.PercentFillPolice = 4; % Percent Lattice is filled with agent type
    p.PercentFillPillars = 1.53; % Percent Lattice is filled with agent type
    p.vision = 4; % Number of spaces an agent can see in each direction
    p.MaxJailTerm = 10; % Max Jail Term in time steps
    p.StartingGovernmentLegitimacy = 0.6; % Starting value (between 0 and 1) for government legitimacy
    p.ChanceFindNVResistor = 40; % Percent chance a police will find a nonviolent resistor who is not visible
    p.ChanceTargetNonviolent = 25; % Percent chance police will target a nonviolent agent
    p.ChanceKillNonviolent = 10; % Percent chance police will attempt to kill nonviolent agent
    p.BackfireCoefficient = .99; % Factor by which government legitimacy decreases
    p.f = .07; % Threshold used for determinging behavior
    p.ProtestCycle = 8.9; % Period of time steps for protesting frequency
    p.ProtestDuration = 1; % Number of timeSteps to protest for each protest cycle
    p.nNV = 1; % Threshold for nonviolent people determining to protest
    p.DefectThreshold = .05; % Threshold for pillars and police to defect
    p.PeerPressureNumber = 3.46; % # of Protestors for PeerPressure to equal .5 (and vision = 7)
    p.DelayStartMax = 1; % Maximum number of time steps for agents to remain inactive
    p.NonviolentSuccessPercent = 10; % Percent of defecting pillars required for nonviolent victory
    p.PillarProxStrategy = 0; % 1 if Activists seek pillars, 0 if Activists do not
    p.DefectThresholdStdDev = .01; %Standard Deviation for DefectThreshold for each Pillar/Policeman
    p.ActivistSearchVision = 10;
    p.ReorderAgentsParam = 1; % Reorders agents, but keeps pillars at the end
    p.PercentFickle = 100; % Percent of Citizens who join and rejoin resistance based on seeing protest. The nonfickle stick remain regardless;
    p.PercentImmediateProtest = 5; % %Percent of Citizens who will immediately protest after joining
end

% Starting Calculations
% Rounding
p.ProtestDuration = round(p.ProtestDuration);
p.ProtestCycle = round(p.ProtestCycle);
p.DelayStartMax = round(p.DelayStartMax);
p.PeerPressureNumber = round(p.PeerPressureNumber);
p.vision = round(p.vision);
p.nNV = round(p.nNV);
% Calculating number of agents for each type
NCitizens = round(p.PercentFillCitizens/100*p.LatticeX*p.LatticeY);
Nactivists = round(p.PercentFillActivists/100*p.LatticeX*p.LatticeY);
NPolice = round(p.PercentFillPolice/100*p.LatticeX*p.LatticeY);
NPillars = round(p.PercentFillPillars/100*p.LatticeX*p.LatticeY);
N = NCitizens + Nactivists + NPolice + NPillars;
% Calculating Variables for Citizens
p.GovernmentLegitimacy = p.StartingGovernmentLegitimacy; % Global government legitimacy value
Ey = exp( a + 0.5 * (b ^ 2) ); 
logy = a+randn(1,NCitizens)*b; % Log of income
y = exp(logy); % Income
hardship = exp(Ey-y)./(1+exp(Ey-y)); % Hardship for Each Citizen
costA = 2 * exp(k2.*p.MaxJailTerm.*y)./(1+exp(k2.*p.MaxJailTerm.*y)) - 1; % Original Cost of being arrested for each citizen, higher for rich
costB = 2 * exp(k2.*p.MaxJailTerm./y)./(1+exp(k2.*p.MaxJailTerm./y)) - 1;% New Cost term, higher for poor
cost = (costA.^4+costB.^4).^.25; % New Cost Term
%cost = costA;

% Initialize Graph and Agents
graph = MakeGraph(p.LatticeX,p.LatticeY,p);
agents.unique = 1:N; %Each Agent has a unique identifying number
[graph, agents] = PlaceAgents(graph,agents); 
agents.type = [repmat({'citizens'},1,NCitizens) repmat({'activists'},1,Nactivists) repmat({'police'},1,NPolice) repmat({'pillars'},1,NPillars)];
agents.alive = ones(1,N); %Every agent starts alive
agents.defect = [nan(1,NCitizens+Nactivists) zeros(1,NPolice+NPillars)];
agents.IndivDefectThreshold = [nan(1,NCitizens+Nactivists) (repmat(p.DefectThreshold,1,NPolice+NPillars)+randn(1,NPolice+NPillars)*p.DefectThresholdStdDev)];
agents.IndivDefectThreshold(agents.IndivDefectThreshold<.01)=.01;
agents.JailTerm = [zeros(1,NCitizens+Nactivists) nan(1,NPolice+NPillars)];
agents.InAction = [zeros(1,NCitizens) nan(1,Nactivists) nan(1,NPolice+NPillars)];
agents.ActiveV = [zeros(1,NCitizens) nan(1,Nactivists+NPolice+NPillars)];
agents.ActiveNV = [zeros(1,NCitizens) nan(1,Nactivists+NPolice+NPillars)];
agents.grievance = [zeros(1,NCitizens) nan(1,Nactivists+NPolice+NPillars)];
agents.hardship = [hardship nan(1,Nactivists+NPolice+NPillars)];
agents.cost = [cost nan(1,Nactivists+NPolice+NPillars)];
agents.MoralPhysicalBarrier = [rand(1,NCitizens)*100 nan(1,Nactivists+NPolice+NPillars)];
agents.InProtest = [zeros(1,NCitizens) zeros(1,Nactivists) nan(1,NPolice+NPillars)];
agents.ProtestCounter = agents.InProtest;
agents.DelayStart = randi(p.DelayStartMax,1,N);
agents.fickle = [rand(1,NCitizens)*100<=p.PercentFickle nan(1,Nactivists+NPolice+NPillars)];
agents.ImmediateProtest = [rand(1,NCitizens)*100<=p.PercentImmediateProtest nan(1,Nactivists+NPolice+NPillars)];
if p.ReorderAgentsParam
    agents = ReorderAgents(agents);
end
AgentHistory(1) = agents; %AgentHistory keeps a record of agent values for each time step

%Initial Plotting
PlotAgents = DeterminePlotAgents(agents); % Sets agent types to be plotted
PlotSymbols = {'g.','c.','cx','bv','kd','ko','k.'}; % Sets symbol colors and shapes to use
PlotSymbolSizes = [20 20 5 5 5 5 5];  % Sets size of symbols to use
PlotAgents2 = DeterminePlotAgents2(agents); % Sets agent types to be plotted
PlotSymbols2 = {'cx','kd','ko','k.'}; % Sets symbol colors and shapes to use
PlotSymbolSizes2 = [5 5 5 5];  % Sets size of symbols to use
if sim.PlotOn
    PlotCustomAgents(graph,PlotAgents,PlotSymbols,PlotSymbolSizes,1); % Makes Graph
    PlotCustomAgents(graph,PlotAgents2,PlotSymbols2,PlotSymbolSizes2,4); hold on % Makes Graph
    figure(2);clf % Clears figure for line plots
    set(gca,'Color',[.8 .8 .8]) % Shades background of line plots to gray
    figure(3);clf % Clears figure for line plots
end

% Run Model
RunCondition = 1; % Ensures initial run of model
ii = 1; % Time Step Variable
victory = 0; % Initializes that there is no victory
MaxInProtestPercentage = 0; % Sets Max Protest Size
ActualMaxInProtestPercentage = 0; % Sets Max Protest Size, This includes people who protest and then get arrested/killed before pillars see them
MaxNVCitizens = 0;
OldProtestPercent = 0;
ProtestPercent = 0;
while(RunCondition)
    % Advance TimeStep
    ii = ii+1;
    
    % Initalizes Protest Count variable
    ProtestCount = 0;
    
    % Run Rules for each Agent
    for jj = 1:length(agents.unique)
        if agents.alive(jj) == 1 && (ii>=(agents.DelayStart(jj)+1)) %Only evaluates for living agents
            % Decrease Protest Counter
            if agents.ProtestCounter(jj) > 0
                agents.ProtestCounter(jj) = agents.ProtestCounter(jj) - 1;
                if agents.ProtestCounter(jj) <= (p.ProtestCycle - p.ProtestDuration)
                    agents.InProtest(jj) = 0;
                end
            end
            % Rule M: Movement Rule, for those not in jail. Pillars don't move
            if ~strcmp(agents.type(jj),'pillars') && ((agents.JailTerm(jj)==0) || isnan(agents.JailTerm(jj)))
                [graph, agents] = RuleM(graph,agents,jj,p);
            end
            % Rule P: Police Arrest or Kill, if they have not defected
            if (agents.defect(jj)==0) && strcmp(agents.type(jj),'police')
                [graph, agents,p] = RuleP(graph,agents,jj,p);
            end
            % Determine Rule Order and run Rules NV and C
            RuleOrder = [1 2] + (agents.ImmediateProtest(jj)==1)*[1 -1];
            for kk = RuleOrder
                if kk == 1
                    % Rule NV: Nonviolent People Protest
                    if ((((agents.JailTerm(jj)==0) && strcmp(agents.type(jj),'activists')) || (agents.ActiveNV(jj) == 1)) && (agents.ProtestCounter(jj)<=0))
                        [graph, agents] = RuleNV(graph,agents,jj,p);
                    end
                    if agents.InProtest(jj)==1  %Adds to Protest Counter if the agent is in protest
                        ProtestCount = ProtestCount + 1;
                    end
                else
                    % Rule C: Determine Citizen will resist violently or
                    % nonviolently, for unjailed citizens
                    if (agents.JailTerm(jj)==0) && strcmp(agents.type(jj),'citizens') && ~(agents.ActiveNV(jj) && ~agents.fickle(jj))
                        [graph, agents] = RuleC(graph,agents,jj,p,k1);
                    end
                end
            end
            % Rule D: Police or Pillars Defect
            if strcmp(agents.type(jj),'police') || strcmp(agents.type(jj),'pillars')
                [graph, agents] = RuleD(graph,agents,jj,p);
            end
            % Decrease JailTerm
            if agents.JailTerm(jj) > 0
                agents.JailTerm(jj) = agents.JailTerm(jj) - 1;
            end

        end
    end
 
    % Pause if requested
    if sim.PauseTime ~= 0
        pause(sim.PauseTime)
    end
    
    % Save state of agents for this time step
    AgentHistory(ii) = agents;
    
    %Plotting
    if sim.PlotOn
        % Create Graph
        PlotAgents = DeterminePlotAgents(agents);
        PlotCustomAgents(graph,PlotAgents,PlotSymbols,PlotSymbolSizes,1);
        title(['Step = ',num2str(ii)])
        % Create Graph of Activists Finding Pillars
        PlotAgents2 = DeterminePlotAgents2(agents);
        PlotCustomAgents(graph,PlotAgents2,PlotSymbols2,PlotSymbolSizes2,4); hold on % Makes Graph
        title(['Step = ',num2str(ii)])
        % Create Line Plots
        x = [ii-1 ii];
        % Determine data to be plotted
        LineData{1} = [sum(strcmp(AgentHistory(ii-1).type,'citizens') & AgentHistory(ii-1).alive & AgentHistory(ii-1).ActiveNV==1) sum(strcmp(AgentHistory(ii).type,'citizens') & AgentHistory(ii).alive & AgentHistory(ii).ActiveNV==1)];
        LineData{2} = [sum(strcmp(AgentHistory(ii-1).type,'citizens') & AgentHistory(ii-1).alive & AgentHistory(ii-1).ActiveV==1) sum(strcmp(AgentHistory(ii).type,'citizens') & AgentHistory(ii).alive & AgentHistory(ii).ActiveV==1)];        
        LineData{3} = [sum(strcmp(AgentHistory(ii-1).type,'activists').*AgentHistory(ii-1).alive) sum(strcmp(AgentHistory(ii).type,'activists').*AgentHistory(ii).alive)];
        LineData{4} = [sum(strcmp(AgentHistory(ii-1).type,'police').*AgentHistory(ii-1).alive) sum(strcmp(AgentHistory(ii).type,'police').*AgentHistory(ii).alive)];
        LineSymbols = {'c-.','r-.','c-','b-'}; %Determin line types
        PlotCustomLines(x,LineData,LineSymbols,2)
        % Create Second Line Plot showing percentage protesting and
        % percentage of pillars defecting
        LineData2{1} = 100*[sum((AgentHistory(ii-1).InProtest==1)&(AgentHistory(ii-1).alive==1))/sum(AgentHistory(ii-1).alive==1) sum((AgentHistory(ii).InProtest==1)&(AgentHistory(ii).alive==1))/sum(AgentHistory(ii).alive==1)]; 
        LineData2{2} = 100*[sum(strcmp(AgentHistory(ii-1).type,'pillars') & AgentHistory(ii-1).alive & AgentHistory(ii-1).defect==1)/sum(strcmp(AgentHistory(ii-1).type,'pillars') & AgentHistory(ii-1).alive) sum(strcmp(AgentHistory(ii).type,'pillars') & AgentHistory(ii).alive & AgentHistory(ii).defect==1)/sum(strcmp(AgentHistory(ii).type,'pillars') & AgentHistory(ii).alive)];
        OldProtestPercent = ProtestPercent;
        LineData2{3} = [OldProtestPercent ProtestCount/sum(agents.alive==1)*100];
        LineSymbols2 = {'b-.','k-','r-.'};
        PlotCustomLines(x,LineData2,LineSymbols2,3)
        xlabel('Time Steps')
        ylabel('%')
        legend('Percent of Population in Protest','Pillar Defection','location','NorthWest')
    end
    
    % Update Max Protest Size
    InProtestPercentage = sum((agents.InProtest==1)&(agents.alive==1))/sum(agents.alive==1)*100;
    if InProtestPercentage > MaxInProtestPercentage
        MaxInProtestPercentage = InProtestPercentage;
    end
    
    % Updates Actual Max Protest Size
    ProtestPercent = ProtestCount/sum(agents.alive==1)*100;
    if ProtestPercent > ActualMaxInProtestPercentage
        ActualMaxInProtestPercentage = ProtestPercent;
    end
        
    % Verify Citizens Resisting in NV
    TotalNVCitizens = sum((agents.ActiveNV==1)&(agents.alive==1));
    if TotalNVCitizens > MaxNVCitizens
        MaxNVCitizens = TotalNVCitizens;
    end
    
    % VictoryCondition
    % Victory set to 1 in previous version of the code that included violent revolution
    if sum(strcmp(agents.type,'pillars').*(agents.defect==1))/sum(strcmp(agents.type,'pillars'))>(p.NonviolentSuccessPercent/100)
        victory = 2; %Sets victory = 2 if nonviolent agents win
    end
    if (sum(strcmp(agents.type,'activists')&(agents.alive==1)) + sum(strcmp(agents.type,'citizens')&(agents.ActiveNV==1)&(agents.alive==1)) <= p.nNV) 
        victory = 3; %Sets victory = 3 if regime wins
    end
    % Run Condition, continue if less than max time steps and no one has won
    RunCondition = ((ii < p.MaxSteps) && (victory == 0));
end

% Output Variables
out.victory = victory;
%out.AgentHistory = AgentHistory;
out.MaxInProtestPercentage = MaxInProtestPercentage;
out.ActualMaxInProtestPercentage = ActualMaxInProtestPercentage;
out.nTimeSteps = ii;
out.MaxNVCitizens = MaxNVCitizens;

if nargin < 2
    % End time measurement
    toc;
end
end

% ***** Supporting Functions *****
% ** Rule Functions **

function [graph, agents] = RuleM(graph,agents,individual,p)
% Movement rule for all agent types
[m, n] = find(graph.agent==individual); % Find agent
%[agents.type(individual) m n]
TinyGraph = CreateTinyGraph(m,n,p.vision,p,graph); % Create Tiny Graph for +/- p.vision around agent

% Find locations with no agents in it
[OpenLocationY, OpenLocationX] = find(TinyGraph.agent==0);
distance = sqrt((OpenLocationX-(p.vision+1)).^2+(OpenLocationY-(p.vision+1)).^2);
OpenNeighborX = OpenLocationX(distance <= p.vision);
OpenNeighborY = OpenLocationY(distance <= p.vision);
if sum(OpenNeighborX)>0 %If there is an open location, move agent
    if ~(strcmp(agents.type(individual),'activists')&&p.PillarProxStrategy~=0)
        ChosenInd = randi(length(OpenNeighborX));
    else
       %Find Pillars
       PillarList = find(strcmp(agents.type,'pillars'));
       if p.PillarProxStrategy==2
            BiggerTinyGraph = CreateTinyGraph(m,n,p.ActivistSearchVision,p,graph);
            [PillarLocationY, PillarLocationX] = find(ismember(BiggerTinyGraph.agent,PillarList));
            NeighborPillars = BiggerTinyGraph.agent(ismember(BiggerTinyGraph.agent,PillarList));    
       else
           [PillarLocationY, PillarLocationX] = find(ismember(TinyGraph.agent,PillarList));
           NeighborPillars = TinyGraph.agent(ismember(TinyGraph.agent,PillarList));    
       end  
       if ~isempty(PillarLocationY)
           if (numel(PillarLocationX)>2) && (p.PillarProxStrategy==2)
               [WeakLinkPillar WeakLinkInd] = min(agents.IndivDefectThreshold(NeighborPillars));
               PillarLocationX = PillarLocationX(WeakLinkInd);PillarLocationY = PillarLocationY(WeakLinkInd);
               PillarLocationX = PillarLocationX - (p.ActivistSearchVision-p.vision);
               PillarLocationY = PillarLocationY - (p.ActivistSearchVision-p.vision);
           end
           OpenNeighborX = [OpenNeighborX; p.vision+1]; %Include Current Location
           OpenNeighborY = [OpenNeighborY; p.vision+1];%Include Current Location
           for jj = 1:length(OpenNeighborX)
               Dist2Pillar(jj,:) = sqrt((PillarLocationX'-OpenNeighborX(jj)).^2+(PillarLocationY'-OpenNeighborY(jj)).^2);
               [MinDist(jj)] = min(Dist2Pillar(jj,:));
           end
           [VeryMinDist VeryMinDistInd] = min(MinDist);
           ChosenInd = VeryMinDistInd;
           if numel(ChosenInd)>1
               ChosenInd = ChosenInd(randi(length(ChosenInd)));
           end
       else % If no Pillars
           ChosenInd = randi(length(OpenNeighborX));
       end
    end
    ChosenSpotX = OpenNeighborX(ChosenInd);
    ChosenSpotY = OpenNeighborY(ChosenInd);
    ChosenM = mod(ChosenSpotY - (p.vision+1)+m,p.LatticeY);ChosenM(ChosenM==0)=p.LatticeY;
    ChosenN = mod(ChosenSpotX - (p.vision+1)+n,p.LatticeX);ChosenN(ChosenN==0)=p.LatticeX;
    graph.agent(ChosenM,ChosenN)=individual; %Place agent in new location in graph
    if ~((ChosenM == m) && (ChosenN == n)) 
        graph.agent(agents.yloc(individual),agents.xloc(individual))=0; %Remove agent from old location in graph
    end
    agents.xloc(individual) = ChosenN; %Update agent variable in X
    agents.yloc(individual) = ChosenM; %Update agent variable in Y
end
end

function [graph,agents,p] = RuleP(graph,agents,individual,p)
% Rule P for policemen, in which they try to arrest or kill
% Determine Neighboring agents and count the relevant agent types
NeighborAgents = DetermineNeighbor(graph,agents,individual,p);
nNVall = sum((strcmp(agents.type(NeighborAgents),'activists').*(agents.JailTerm(NeighborAgents)==0))|(strcmp(agents.type(NeighborAgents),'citizens').*(agents.ActiveNV(NeighborAgents)==1)));
nNV = sum(agents.InProtest(NeighborAgents) == 1); 

FindResistor = rand*100;
CanFindNVResistor = FindResistor < p.ChanceFindNVResistor; % True if police can find hidden nonviolent resistor
TargetNonViolent = rand*100;
if (nNVall > 0) && (TargetNonViolent < p.ChanceTargetNonviolent) && (CanFindNVResistor || nNV>0) % Determine if they target a nonviolent person
    % Targets nonviolent agent
    if (nNV>0) 
        NeighborNonViolent = NeighborAgents(agents.InProtest(NeighborAgents)==1); % Only Target Protestors
    else
        NeighborNonViolent = NeighborAgents(((strcmp(agents.type(NeighborAgents),'activists').*(agents.JailTerm(NeighborAgents)==0))|(strcmp(agents.type(NeighborAgents),'citizens').*(agents.ActiveNV(NeighborAgents)==1))));
    end
    NonViolentTarget = NeighborNonViolent(randi(length(NeighborNonViolent)));
    KillNonViolent = rand*100;
    if KillNonViolent < p.ChanceKillNonviolent % Determine if they try to arrest or kill the nonviolent person
        [graph, agents] = Kill(graph,agents,NonViolentTarget); % Attempts to kill a nonviolent person are always succesful
        p.GovernmentLegitimacy = p.GovernmentLegitimacy*p.BackfireCoefficient; 
    else
        agents = Arrest(agents,NonViolentTarget,p.MaxJailTerm); %Arrests agent
    end
end
end

function [graph, agents] = RuleC(graph,agents,individual,p,k1)
% Rule C for citizens to determine if they will be actively resisting and
% if they choose violence or nonviolence
% Calculates individual agent grievance
agents.grievance(individual) = agents.hardship(individual)*(1-p.GovernmentLegitimacy);
% Determine Neighboring agents and count the relevant agent types
NeighborAgents = DetermineNeighbor(graph,agents,individual,p);
nV = sum(agents.InAction(NeighborAgents) == 1);
nPo = sum(strcmp(agents.type(NeighborAgents),'police'));
nAct = sum(strcmp(agents.type(NeighborAgents),'activists').*(agents.JailTerm(NeighborAgents)==0));
nNVC = sum(strcmp(agents.type(NeighborAgents),'citizens').*(agents.ActiveNV(NeighborAgents)==1));
nInProtest = sum(agents.InProtest(NeighborAgents) == 1);
% *** Very Important **** Choose appropriate ArrestProbability
%ArrestProbability = (1 - exp(-k1*floor(nPo/(1+nNVC+nV+nAct)))); % Arrest Probabilty in NetLogo of Moro ABM
ArrestProbability = (1 - exp(-k1*(nPo/(1+nNVC+nV+nAct)))); % Arrest Probabilty in report of Moro ABM
%NVPeerPressure = .5/p.PeerPressureNumber*(nAct+nNVC); NVPeerPressure(NVPeerPressure>1)=1;
NVPeerPressure = .5/p.PeerPressureNumber*(nInProtest); NVPeerPressure(NVPeerPressure>1)=1;
NVChoice = agents.grievance(individual)*NVPeerPressure - ArrestProbability*agents.cost(individual);

% Agent decides to be part of nonviolent resistance or not
if NVChoice > p.f  % Threshold to become nonviolent
    agents.ActiveNV(individual) = 1;
else
    agents.ActiveNV(individual) = 0;
end
end

function [graph, agents] = RuleD(graph,agents,individual,p)
% Rule D for policeman and pillars to defect
% Determine Neighboring agents and count the relevant agent types
NeighborAgents = DetermineNeighbor(graph,agents,individual,p);
nAct = sum(strcmp(agents.type(NeighborAgents),'activists').*(agents.InProtest(NeighborAgents)==1));
nNVC = sum(strcmp(agents.type(NeighborAgents),'citizens').*(agents.InProtest(NeighborAgents)==1));
NVratio = (nNVC + nAct) / length(NeighborAgents);
if (NVratio > agents.IndivDefectThreshold(individual))  % Defects based on DefectThreshold
    agents.defect(individual) = 1; 
end
end

function [graph, agents] = RuleNV(graph,agents,individual,p)
% Rule NV for nonviolent agents to protest
% Determine Neighboring agents and count the relevant agent types
NeighborAgents = DetermineNeighbor(graph,agents,individual,p);
PoV = sum(strcmp(agents.type(NeighborAgents),'police'));
NVCV = sum(agents.ActiveNV(NeighborAgents) == 1);
AV = sum(strcmp(agents.type(NeighborAgents),'activists') & (agents.JailTerm(NeighborAgents)==0));
% Select rule
%if (PoV == 0) || (AV + NVCV)/PoV > p.nNV % Like InAction Rule, ratio of resistors to police
%if (AV + NVCV)/length(NeighborAgents) > p.nNV % ratio of resistors to all people
if (AV + NVCV) >= p.nNV % ratio of resistors to all people
    agents.InProtest(individual) = 1;
    agents.ProtestCounter(individual) = p.ProtestCycle;
end
end

% ** Supporting Functions to Rules **

function [graph, agents] = Kill(graph,agents,Victim)
% Kills an agent
agents.alive(Victim) = 0; % Alive variable set to 0
graph.agent(agents.yloc(Victim),agents.xloc(Victim)) = 0; % Removed from graph
agents.xloc(Victim) = nan; %Location data removed from agent variable
agents.yloc(Victim) = nan; %Location data removed from agent variable
end

function agents = Arrest(agents,ViolentTarget,MaxJailTerm)
% Arrests an agent
agents.JailTerm(ViolentTarget) = randi(MaxJailTerm); % Jail term assigned
agents.InAction(ViolentTarget) = 0; % Agent cannot attempt to kill
agents.InProtest(ViolentTarget) = 0; % Agent cannot protest
if strcmp(agents.type(ViolentTarget),'citizens') % If a citizen, then they no longer are actively resisting
    agents.ActiveNV(ViolentTarget) = 0;
end
end

function NeighborAgents = DetermineNeighbor(graph,agents,individual,p)
% For a given agent, i.e. individual, this function determines the
% neighboring agents in the radius specified by p.vision
vision = p.vision;
[m, n] = find(graph.agent==individual); 
TinyGraph = CreateTinyGraph(m,n,p.vision,p,graph);
[LocationY LocationX] = find(TinyGraph.agent~=0 & TinyGraph.agent~=individual);
distance = sqrt((LocationX-(p.vision+1)).^2+(LocationY-(p.vision+1)).^2);
NeighborX = LocationX(distance <= p.vision);
NeighborY = LocationY(distance <= p.vision);
NeigbhorAgentLoc = (NeighborX-1)*(p.vision*2+1)+NeighborY;
NeighborAgents = TinyGraph.agent(NeigbhorAgentLoc);
end

% ** Graph Functions **

function [graph] = MakeGraph(LatticeX,LatticeY,p)
% Make Graph
graph.x = repmat(1:LatticeX,LatticeY,1);
graph.y = repmat((1:LatticeY)',1,LatticeX);
graph.agent = zeros(LatticeY,LatticeX);
% Make UberGraph, i.e. graph extended by vision variable in each direction
% [ydim, xdim]  = size(graph.x);
% graph.ubergraph.x = [ repmat(-p.vision+1:xdim+p.vision,p.vision,1)
%     repmat(-p.vision+1:0,ydim,1) graph.x repmat(xdim+1:xdim+p.vision,ydim,1)
%     repmat(-p.vision+1:xdim+p.vision,p.vision,1)];
% graph.ubergraph.y = [ repmat((-p.vision+1:ydim+p.vision)',1,p.vision)...
%     [repmat((-p.vision+1:0)',1,xdim); graph.y; repmat((ydim+1:ydim+p.vision)',1,xdim)]...
%     repmat((-p.vision+1:ydim+p.vision)',1,p.vision,1)];
% graph.ubergraph.agent = MakeTorusMatrix(graph.agent,p);
% graph.ubergraph.xprime = MakeTorusMatrix(graph.x,p);
% graph.ubergraph.yprime = MakeTorusMatrix(graph.y,p);
end

function [graph, agents] = PlaceAgents(graph,agents)
% Creates initial placement of the agents in the graph
for ii = 1:max(agents.unique)
    OpenLocation = find(graph.agent(:)==0);
    AgentLocation = OpenLocation(randi(length(OpenLocation)));
    graph.agent(AgentLocation) = ii; %Updates graph variable
    agents.xloc(ii) = graph.x(AgentLocation); %Updates agent variable
    agents.yloc(ii) = graph.y(AgentLocation); %Updates agent variable
end
end

function [out] = MakeTorusMatrix(in,p)
% Increases the size of a graph matrix, extending it in all directions by vision variable
temp = [in(:,end-p.vision+1 : end) in in(:,1:p.vision)];
out = [ temp(end-p.vision+1:end,:); temp; temp(1:p.vision,:)]; 
end

function TinyGraph = CreateTinyGraph(m,n,vision,p,graph)
% Creates Graph +/- p.vision around selected agent,
TinyGraph.x = repmat(mod(n-vision : n+vision,p.LatticeX),2*vision+1,1);
TinyGraph.x(TinyGraph.x==0)=p.LatticeX;
TinyGraph.y = repmat(mod(m-vision : m+vision,p.LatticeY)',1,2*vision+1);
TinyGraph.y(TinyGraph.y==0)=p.LatticeY;
tempind = (TinyGraph.x-1)*p.LatticeY+TinyGraph.y;
TinyGraph.agent = graph.agent(tempind);
end

% ** Plotting Functions **

function PlotCustomAgents(graph,PlotAgents,PlotSymbols,PlotSymbolSizes,PlotHandle)
% Creates a plot of the graph
figure(PlotHandle)
clf
set(gca,'Color',[.8 .8 .8]); %Sets background to gray
hold on
for ii = 1:length(PlotSymbols) %For loop for plotting each agent as specified in PlotAgents
    inds = ismember(graph.agent,PlotAgents{ii});
    plot(graph.x(inds),graph.y(inds),PlotSymbols{ii},'MarkerSize',PlotSymbolSizes(ii))
end
hold off
xlim([0 max(graph.x(:))+1])
ylim([0 max(graph.y(:))+1])

end

function PlotAgents = DeterminePlotAgents(agents)
% Determines the agents to be plotted in the graph
PlotAgents{1} = agents.unique(strcmp(agents.type,'citizens')); %Citizens
PlotAgents{2} = agents.unique(strcmp(agents.type,'citizens')&(agents.ActiveNV==1)); %Nonviolent Citizens
PlotAgents{3} = agents.unique(strcmp(agents.type,'activists')); %Activists
PlotAgents{4} = agents.unique(strcmp(agents.type,'police')); %Police
PlotAgents{5} = agents.unique(strcmp(agents.type,'pillars')); %Pillars
PlotAgents{6} = agents.unique(agents.InProtest==1);
PlotAgents{7} = agents.unique(agents.JailTerm > 0);
end

function PlotAgents = DeterminePlotAgents2(agents)
% Determines the agents to be plotted in the graph
%PlotAgents{1} = agents.unique(strcmp(agents.type,'citizens')); %Citizens
%PlotAgents{2} = agents.unique(strcmp(agents.type,'citizens')&(agents.ActiveNV==1)); %Nonviolent Citizens
PlotAgents{1} = agents.unique(strcmp(agents.type,'activists')); %Activists
%PlotAgents{6} = agents.unique(strcmp(agents.type,'police')); %Police
PlotAgents{2} = agents.unique(strcmp(agents.type,'pillars')); %Pillars
PlotAgents{3} = agents.unique(agents.InProtest==1);
PlotAgents{4} = agents.unique(agents.JailTerm > 0);
end

function PlotCustomLines(x,LineData,LineSymbols,PlotHandle)
% Creates a line plot agent counts at each time step
figure(PlotHandle)
hold on;
PlotString = [];
for ii = 1:length(LineData)
    PlotString = [PlotString,'[',num2str(x),'], LineData{',num2str(ii),'}, LineSymbols{',num2str(ii),'},'];
end
PlotString = PlotString(1:end-1);
eval(['plot(',PlotString,')'])
end

% ** Other Functions **

function out = ReorderAgents(in)
% Function reorders agents so the actions of different agent types are
% interspersed.
% Determine non-pillar agents
NonPillarAgents = ~strcmp(in.type,'pillars');
% Create New Order
MaxNonPillarAgent = max(NonPillarAgents.*in.unique);
NewOrder = [randperm(MaxNonPillarAgent) [MaxNonPillarAgent+1:max(in.unique)]];
% Prepare Loop for creating new agent structure
AgentFields = fieldnames(in);
out = [];
ExcludeFields = {'unique','xloc','yloc'};
for ii = 1:length(AgentFields);
    ThisField = AgentFields{ii};
    OldField = getfield(in,ThisField);
    if sum(strcmp(ThisField,ExcludeFields))==0
        NewField = OldField(NewOrder); %Reorder structure fields
    else
        NewField = OldField; %Do Not Reorder excluded fields as that would mess up graph
    end
    out = setfield(out,ThisField,NewField);
end
end
